arXiv:1504.07995v2 [astro-ph.EP] 12 Oct 2015 


Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 14 October 2015 (MN style file v2.2) 


An Empirically Derived Three-Dimensional Laplace 
Resonance in the Gliese 876 Planetary System 

Benjamin E. Nelson^’^, Paul M. Robertson^’^, Matthew J. Payne^, Seth M. Pritchard^, 
Katherine M. Deck^, Eric B. Eord^’^, Jason T. Wright^’^, Howard T. Isaacson® 

^ Center for Exoplanets and Habitable Worlds, The Pennsylvania State University, 525 Davey Laboratory, University Park, PA, 16802, USA 
^Department of Astronomy & Astrophysics, The Pennsylvania State University, 525 Davey Laboratory, University Park, PA 16802, USA 
^ Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA 

'^Department of Physics & Astronomy, University of Texas San Antonio, UTSA Circle, San Antonio, TX 78249, USA 
^Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, CA 91101, USA 
^Department of Astronomy, University of California, Berkeley, Berkeley, California 94720, USA 


14 October 2015 


ABSTRACT 

We report constraints on the three-dimensional orbital architecture for all four planets 
known to orbit the nearby M dwarf Gliese 876 based solely on Doppler measurements 
and demanding long-term orbital stability. Our dataset incorporates publicly available 
radial velocities taken with the ELODIE and CORALIE spectrographs, HARPS, and 
Keck HIRES as well as previously unpublished HIRES velocities. We first quantita¬ 
tively assess the validity of the planets thought to orbit GJ 876 by computing the 
Bayes factors for a variety of different coplanar models using an importance sampling 
algorithm. We find that a four-planet model is preferred over a three-planet model. 
Next, we apply a Newtonian MCMC algorithm to perform a Bayesian analysis of the 
planet masses and orbits using an n-body model in three-dimensional space. Based on 
the radial velocities alone, we find that a 99% credible interval provides upper limits 
on the mutual inclinations for the three resonant planets ($cb < 6.20° for the c and b 
pair and <I>be < 28.5° for the b and e pair). Subsequent dynamical integrations of our 
posterior sample find that the GJ 876 planets must be roughly coplanar ($cb < 2.60° 
and $be < 7.87°), indicating the amount of planet-planet scattering in the system 
has been low. We investigate the distribution of the respective resonant arguments 
of each planet pair and find that at least one argument for each planet pair and the 
Laplace argument librate. The libration amplitudes in our three-dimensional orbital 
model supports the idea of the outer-three planets having undergone significant past 
disk migration. 

Key words: planets and satellites: dynamical evolution and stability - planets and 
satellites: individual - planets and satellites: formation - techniques: radial velocities 
- methods: numerical - methods: statistical 


1 INTRODUCTION 

Gliese 876 (=G J 876) is a 0.37 Mq M4V star 
llvon Braun et all l2014h hosting four known planets. This 
remarkable system represents several milestones: the first 
detection of a planet around an M-d warf (GJ 876 b) 
llMarcv et al.l 1 19981 : iDelfosse et al.l Il998l ). the first exam- 
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pie of multi-plan et system orbiting in a mean-motion res¬ 
onance (MMR) llMarcv et al.l l200lh . the first exar nple o f 
an MMR chain amongst three planets llRivera et al .1120101 ) . 
and the closest mu lti-planet exosystem to date (4.689 pc, 
Ivan LeeuwerJ (l2007l )). 

The star has a lengthy Doppler (or radial velocity, RV) 
history spanning two decades and multiple ob serving sites. 
Plane t b was detected contemporaneously by iMarcv et al.l 
lll998l) using the Lick Hamilton Spectrograph and Keck 
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HIRES and IPelfosse et alj (Il998l 'l using the ELODIE and 
CORALIE spectrographs. Both estimated a moderately ec¬ 
centricity for b (~0.3) and an orbital period of 61 days 
for this gas giant fro m their radia l velo city model. With 
more RV observations, iMarcv et al.l (1200 ll i uncovered a sec¬ 
ond gas giant, c, orbiting near 30 days. This planet’s RV 
signature previously masqueraded as a larger eccentricity 
for planet b due t o the near 2:1 period commensurabil- 
ity o f their orbits (lAnglada-Escude et al.ll2010l: lEord et al.l 
M). As the Keck RV dataset grew. iRivera et al.l (l2005l l 
revealed a third planet d orbiting near 1.9 days and 
was the lowest mass exoplanet around a main-sequence 
star at the time (m sini=5.89M@). Ph otometric measure¬ 
ments showed plane t d did not transi t (|Rivera_eLal, II 2 OO 5 I: 
Laughlin et al.l[2005l: IShankland et al.|[2CI06l : iKammer et al.l 


20lj i. Correia et al.l ( 2010l ') published new HARPS RVs 


which by themselves could constrain the mutual inclina- 
tion between plan ets c and b. Around the same time, 
iRivera et al. liioIO) published new Keck RVs which showed 
an additional RV signal around 124 days, dubbed planet e. 
Numerically integrating their solutions beyond the last ob¬ 
servation, the outer three planets (c, b, and e) appear to be 
in a Laplace resonance, much like the three closest Galilean 
moons orbiting Jupiter. Other studies have placed limits 
on the existence of additional planets and massive com¬ 
panions in the system through observations dLeinert etliD 
I 1997 I : iLuhman fc Javawardhan3l2002l : IPatience et al.ll2002l l 
and considerations of long-term dynamical stability 


dJones et al.l 2001: Ji & Liul 20061: 

Rivera & Haghighipouil 

I 2 OO 7 I: iGerlach & Haghighipouil 2012 

)• 


For RV systems, we only observe the component of the 
planetary induced stellar wobble projected onto our line- 
of-sight. Most of the time, there is a degeneracy between 
the true mass (m) and on-sky inclination (i), where an 
edge-on system is isys = 90°, so we can only place a lower 
limit on the orbiting companion’s mass. However, if the self¬ 
interactions in a multiple planet system are strong, the RV 
model becomes sensitive to the true masses of the planets, 
thereby breaking the m sin isys degeneracy. There are many 
RV systems where the true masses can be me aningfully con¬ 
strain ed, including HD 200964, 24 Sextantis J Johnson et al.l 
l201lh . HD 829 43 (iTan et al]l2013h . and other dynamically 
active systems (jVeras fc FoiHjlgOljjL _ 

For GJ 876 c and b, iLaughlin fc Ghamberj ll200lll and 
IRivera fc Lissau^ ll200ll ~) performed self-consistent Newto¬ 
nian fits and constrai n the planetar y mass es and their on-sky 
coplanar inclination. IRivera et al.l (l2005l l found an on-sky 
inclination for the t hree-planet system (isvs = 50° ±3°), 
assuming coplanarity. Bean_fc SeifahrtI (|2009l 'l combined the 
RVs from IRivera et a hl (|2005l ) and Hub ble Space Telescope 
astrometry from Benedict et al.l 1 I 2 OO 2 II and measured the 
mutual inclination between c and b, where 

cos <E>cb = cos ic cos ib J- sin ic sin ib cos (He — ^b) , (1) 

H is the longitude of ascending node, and c and b de¬ 
note the two planets considered. They find iheb = 5.0± ?’^o. 
Based on the HARPS RVs alone, ICorreia et al.1 1 I 2 OIOII find 
a lower value (4?cb = 1.00°). With their four-planet model, 
IRivera et al.l (l201(tl finds a best fit value 4>cb = 3.7°. All of 
the above values are generally consistent with a nearly copla¬ 


nar s ystem. When incorporating the effects of correlated 
noise, iBaluevI (l201ll 'l places an upper limit of "heb = 5 — 15°. 

Arguably the most studied exoplanet system display¬ 
ing a MMR, the GJ 876 system has been a prime testbed 
for dynamical and planet formation theory for the past 
decade. The 2:1 MMR of the c and b pair is of par¬ 
ticular interest, since the libration amplitude of the res¬ 
onant arguments are a valuable ind icator of the sys¬ 
tem’s long-term d y namical stability jKinoghii^^^^^^ 

2001; [^^t^2j2002jGoidzie^kre^M2|^^3|^^^^^^^3 


2003; Haghigliipom et al.1 2003l :~ Zhou fc Sunl ~ ^03 L The 

most likely mechanism for the planets’ current orbital 
periods and eccentricities is through a combination of 
planet-planet interactions and disk migration, while most 
studies have focused o n migration throug h a gas disk 
llSnellgro ve et al.l I2001I: iMurrav et al.1 12002: iL ee fc Pe^ 
2OO2I: iKlev et al.l l2004 iLed 120*^ iKlev et al.l 12005 


Thomm esI I2OO5 I: Thommes et al.l l2( 0^;_jLee_j£^rhonin^ 

200d; IPodlewska-Gaca et al.1 I2OI2I : Batygin fc Morbidellil 
2OI3I : iLega et al.1 20131 '). fSemil-analvtic models describing 
the s ecular evolution of the system have also been devel¬ 
oped (iBeauge et al.|[200^ : IVerasll2007l 'l. 

Other studies of the GJ 876 system include the search 
for debris disks, thermal properties of d, and interior struc¬ 
ture models of d (IRivera et al.ll20l^ . and references therein). 
More recent analyses include how the star’s UV radi - 
ation field affects habitability (iFrance et al.l I 2 OI 2 I . l2013f l 
and the detect ability of additional hypothesi zed planets 
in the system dGerlach fc Haghighipouil I 2 OI 2 IL Although 
the system’s orbital architecture appears atypical, GJ 876 
does fit into a larger portrait of exoplanet systems around 
M dwarfs, includin g population statistics of giant planets 
dMontet et al.ll2 014l~l and the co planaritv of multi-planet sys¬ 
tems ( Ballard fc Johnsordl201 j l. 

RV surveys have discovered a couple dozen strongly- 
interacting multi-planet systems, motivating the need for 
analysis pro cedures to inc o rporat e a self-consistent Newto¬ 
nian model. iNelson et al.l d2014al ) developed a Newtonian 
MCMG algorithm which has been successful in analyzing 
systems that re quire a large numbe r of model parameters 
(e.g. 55 Gancri: [kelson et al.ll2014blL The GJ 876 observa¬ 
tions have a similar observing baseline and require almost 
as many model parameters, making it compatible with such 
an algorithm. 

In this work, we present a detailed characterization of 
the orbits and masses of the GJ 876 planets employing a full 
three-dimensional orbital model used to fit the Doppler ob¬ 
servations. In we describe the RV observations made with 
multiple spectrographs (ELODIE, GORALIE, HARPS, and 
HIRES), including a new set of HIRES measurements. In 
Sj3l we describe our orbital and observational model and in¬ 
vestigate the effects of correlated noise in the observations. 
In Sj4l we report our results for a coplanar orbital model. 
Before advancing to a more complex orbital model, we eval¬ 
uate the evidence for planets d and e by computing Bayes 
factors as described in ()3.4I In |j5l we report the results of 
the fitting and n-body simulations for a three-dimensional 
orbital model. In we investigate the evolution of the res¬ 
onant angles and statistics regarding their libration ampli¬ 
tude. In Appendix m we test for possible observational bi- 
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ases in these results. We conclude with a discussion of the 
key results and the applications of our posterior samples in 

El 


2 OBSERVATIONS 

Our dataset includes publicly-available RVs from four differ¬ 
ent instruments. We include 46 ELODIE, 40 CORALIE, and 
52 HARPS observations (jCorreia et al.l[20inl l. The inclusion 
of these data extends our observing baseline to roughly 580 
days before the first Keck HIRES observation. The 162 Keck 
observ ations are reduced by the Carnegie Planet Search 
group ll Rivera et al. IH). These will be referred to as the 
Carnegie RVs henceforth. 

Our analysis also includes 67 additional Doppler mea¬ 
surements from HIRES reduced by the California Planet 
Search group (Table [T]). These will be referred to as the 
California RVs henceforth. We measured r elative RVs of GJ 
876 with the HIRES echelle spectrometer llVogt et al.lll99^ i 
on the 10-m Keck I telescope using standard procedures. 
Most observations were made with the B5 decker (3.5 x 
0.86 arcseconds). Light from the telescope passed through a 
glass cell of molecular iodine cell heated to 50° C. The dense 
set of molecular absorption lines imprinted on the stellar 
spectra between 5000-6200 A provide a robust wavelength 
scale against which Doppler shifts are measured, as well as 
strong constraints on the instrumenta l profi l e at the time 
of ea ch observation (iMarcv fc ButleJ Il992l : IValenti et al.l 
Il995h . We also obtained five iodine-free “template” spectra 
using the B1 decker (3.5 x 0.57 arcseconds). These spectra 
were de-convolved using the instrumental profile measured 
from spectra of rapidly rotating B stars observed immedi¬ 
ately before and after through the iodine cell. We measured 
high-precision relative RVs using a forward model where the 
de-convolved stellar spectrum is Doppler shifted, multiplied 
by the normalized high-resolution iodine transmission spec¬ 
trum, convolved with an instrumental profile, and matched 
to the observed spectra using a Levenb erg-Marquardt algo - 
rithm that minimizes the statistic jButler et al.|[l99^ . 
In this algorithm, the RV is varied (along with nuisance pa¬ 
rameters describing the wavelength scale and instrumental 
profile) until the minimum is reached. Each RV uncer¬ 
tainty is the weighted error on the mean RV of ~700 spec¬ 
tral chunks that are separately Doppler analyzed. These un¬ 
certainty estimates do not account for potential systematic 
Doppler shifts from instrumental or stellar effects. 


3 METHODS 

We characterize the masses and orbits of the GJ 876 plan¬ 
ets using the GPU version of the Radial velocity Using 
N-body Differ ential evolution Mar kov Chain Monte Carlo 
code, RUNDMC ll Nelson et aLll^ldal ). which incorporates the 
Swarm-NG fram ework to integrate planetary systems on 
graphics cards llPindar et al.l l2013l ) . RUN DMC has analyzed 
the 55 Cancri planetary system, which required a high¬ 
dimensional (^40) model. Compared with 55 Cancri, the 
problem of GJ 876 seems to be similarly challenging but 


more computationally tractable due to having fewer obser¬ 
vations, a shorter observing baseline, a larger integration 
timestep, and a lower dimensional model (e.g. one less planet 
to account for). On the other hand, the presence of extremely 
strong planet-planet interactions can result in a challenging 
posterior distribution that could be more difficult to sample 

from. _ _ 

Considering the lessons from I Nelson et al.l ll2014al i re¬ 
garding how to explore parameter space efficiently, we set 
the following algorithmic parameters for RUNDMC: richains = 
300, a-y = 0.05, and MassScaleFactor=1.0. To accommo¬ 
date the inner-most planet’s 1.9 day orbital period, we set 
our integration timestep to roughly 17 minutes and use th e 
time-symmetrized Hermite integrator (iKokubo et al.l[l998l L 
We begin by apply ing RUNDMC to the Carnegie RVs 
from iRivera et akl (l2010li for a four-planet model. We sam¬ 
pled from the Markov chains after they had burned-in suffi¬ 
ciently. In RUNDMC, we specify the orbital parameters at the 
epoch of the first observation in our dataset. For the full 
RV dataset, this happens to be an ELODIE observation. To 
generate the initial states of Markov chains for analyzing the 
full dataset, we start from a posterior sample based on the 
Carnegie RVs and rewind each planet’s argument of perias- 
tron (oj) and mean anomaly (M) according to the preces¬ 
sion rate (th), orbital periods, and time diff erence between 
the fir st ELODIE and HIRES observations. iLaughlin et al.l 
1I2OO5II find the joint line of apses for the c and b pair to 
be precessing at an average rate oi vj — —41° yr“^. After 
burning-in, we randomly sample from our Markov chains to 
use as our set of initial conditions for long-ter m stability 
simu lations using the MERCURY hybrid integrator llChamberj 
Il999lf . The solutions vetted for stability are used as initial 
conditions for more restricted RUNDMC runs to be explained 
in El 

3.1 Model Parameters 

We characteriz e the system model wit h a fixed star mass 
(M* =0.37Mr;>: lvon Braun et al.l (l2014ll i. plus each planet’s 
mass (m), semi-major axis (o), eccentricity (e), inclination 
(i), argument of periastron (w), longitude of ascending node 
(D), and mean anomaly (M) at our chosen epoch (first 
ELODIE observation) for each planet, plus the RV zero- 
point offsets (C) and jitters (i.e. unmodeled instrumental 
and astrophysical noise, Oju) for each observatory. We re¬ 
port the orbital periods (P) based on Kepler’s Third Law 
and each body’s m and a based in a Jacobi coordinate 
system. Planet masses and semi-major axes can be readily 
rescaled to account for any updates to the stellar mass. 

The GJ 876 planets are well approximated by a coplanar 
system, i.e. D = 0° and i is the same for all planets. Due 
to the symmetrical nature of a radial velocity system on the 
sky, a planetary system at isys is indistinguishable from one 
at an inclination of 180° — isys- So in El we restrict the 
planets to a coplanar system that can take on any value of 
0° < isys < 90°. Similarly, the entire system can be rotated 
about the line-of-sight. In El the individual i and Q. for each 
planet become free parameters in our model that are fit to 
the observations. We set Q,d = 0° to ground our coordinate 
system. 
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Table 1. New Keck HIRES velocities for GJ 876. 


BJD-2450000. [days] 

Radial Velocity [ms ^ ] 

Uncertainty [ms ^ ] 

3301.808032 

-33.89 

2.06 

3301.816875 

-36.49 

2.37 

3301.822986 

-28.19 

2.60 

3301.871412 

-39.54 

1.99 

3302.722523 

-91.03 

1.91 

3302.728866 

-89.20 

2.02 

3302.735567 

-81.93 

1.92 


Table [T] is presented in its entirety as Supporting Information with the online version of the article. This stub table is shown for 

guidance regarding its form and content 


3.2 Model of Observations 

RUNDMC allows for fitting multiple zero-point offsets and 
magnitudes of jitter (e.g. combined astrophysical and/or in¬ 
strumental noise). HIRES received a CCD upgrade and new 
Doppler reduction process in August 2004 (JD 2453241.5). 
The Carnegie time series was split based on this pre- and 
post-upgrade era, which was modeled by two zero-point off¬ 
sets, but the entire dataset is still modeled with one jitter 
term. The California RVs also have a separate offset and one 
jitter term. The other instruments (ELODIE, CORALIE, 
HARPS) were each modeled by one RV zero-point offset 
and jitter term. In total, we account for six offsets and five 
jitter parameters. 

Our datase t does not include the HST astrometry from 
[Benedict et al.l (|2002| ~) or any informative priors regarding 
the inclinations of the GJ 876 planets. 


3.3 Magnitude and Timescale for Correlated 
Noise 

Our likelihood funct ion described in Equations 4 and 5 of 
I Nelson et all (|2014ah assumes that the observational errors 
are uncorrelated, which may not be a sufficient approxi¬ 
mation for high precision RV me a surem ents of stars with 
significant stellar activity. iBaluevI (|201ll ~l showed that cor¬ 
related (=red) noise in the RV data could lead to bias or 
misestimated uncertainty in some of the orbital parame¬ 
ters (e.g., Bd). Recently, some discoveries of planets orbit¬ 
ing M dwarfs have been shown to b e more likely mere arti¬ 
facts resulting from stellar act ivity dRobertson et al]l2014l : 
[Robertson fc Mahadevanll2014h . Therefore, we take a closer 
look at the role of stellar activity in contributing astrophys¬ 
ical noise to RV observations of GJ 876. In 113.41 we perform 
a complementary analysis by computing Bayes factors for a 

finite set of models. _ 

_ High-precision photometry of GJ 876 dRivera et al.[ 

[2005l i revealed a rotation period of 97 ± 1 days. Examin¬ 
ing the variability of the activity-sensitive Ha and Na I D 
absorption lines in the publicly-available HARPS spectra of 
the staiQ, we confirm the rotation period, finding an Ha 


^ Based on data obtained from the ESO Science Archive Facility 
under request number 129169 



Figure 1. (Top) Genera lized Lomb-Scargle periodogram 
dZechmeister &: Kurster|[2009[i of the Ha stellar activity index for 
GJ 876, as measured from the publicly available HARPS spectra. 
We mark the peak at P = 95 days—which we adopt as the stellar 
rotation period-and its 1-year alias at 128 days. (Bottom) Resid¬ 
ual periodogram of Ha after modeling and removing a sinusoid at 
the rotation period. The dashed lines indicate the power required 
for a false alarm probability of 1 p ercent according to Equation 
24 of [Zechmeister fc Kiirsten ll200d ~). 


periodicity of 95±1 days (Figure[T] top panel). The appear¬ 
ance of the rotation period in the photometry suggests the 
pres ence of starspots, which can affect the m easured RVs 
fe.g. [Boisse et all 2011) : Dumusaue et al.l[2014[) . The results 
from BaluevI l[201lj) . Rivera et al.l 1 20051 ). and our own ex¬ 
amination of the spectral activity tracers suggest a complete 
treatment of activity will yield only marginal improvements 
in the accuracy of the RV model, and is not necessary to 
achi eve the goals of thi s study. 

[Boisse et aL[ |[201lj) showed that rotating starspots will 
induce RV periodicities at the stellar rotation period Piot 
and its integer fractions (Prot/2, Prot/3, etc.). Given that 
none of the planets in the GJ 876 system have periods near 
the rotation period or its integer fractions, we conclude that 
rotating starspots have not produced large-amplitude peri¬ 
odic RV signals such as might be misinterpreted as planet 
candidates. Although the RV amplitude of the outermost 
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planet e (3.5ms“^) is closest to the amplitudes expected for 
activity signals produced by a chromospherically quiet star 
such as GJ 876, we are unaware of any physical mechanism 
that would create an RV signal on this timescale. 

Our computed stellar rotation period of 95 ± 1 days 
beating with one Earth year produces an alias of around 
128 days, which is worryingly close to the orbital period 
of planet e (~124 days). However, subtracting the rotation 
signal causes the 128-day peak to disappear as well, leading 
us to suspect that this signal is an alias (Figure [U bottom 
panel). If the 124-day signal is an 1-year alias of the rotation 
perio d, then we would als o see the 95 d ay sig nal in the RV 
data. I Rivera et al. I ll201(]ll and ISaluevI (1201 il l did not find 
the stellar rotation signal in the RV datasets they analyzed. 


3.4 Computing Bayes Factors for Model Selection 


When Doppler observations of a system with multiple 
strongly-interacting planets, the shape of the posterior dis¬ 
tribution is often challenging to sample from efficiently. A 
general Bayesian approach of performing model comparison 
is to compute the fully marginalized likelihood, sometimes 
called the evidence, for each model. Formally, the evidence 
is the probability of generating the observed radial veloc¬ 
ity dataset d assuming some underlying model M that is 
parameterized by 6, 

p{d\M) = J pid\e,M)pie\M)de ( 2 ) 

where p{d\6, M) is the likelihood function and p{9\M) is the 
prior probability distribution. While the value of p{d\Ai) 
is not useful by itself, the ratio of two evidences for two 
competing models Mi and M 2 , yields the Bayes factor 

BF^p(j-M2) (3^ 

p{d\Mi) 

that provides a quantitative measure of which model is pre¬ 
ferred and to what degree. 

Marginal likelihoods are notoriously difficult to com¬ 
pute. Integrating Equation [2] analytically may be possible 
for some idealized problems with one to a few dimensions, 
but the model required to describe a 3-1- planet system needs 
roughly 20 or more parameters. Numerical integration tech¬ 
niques such as Monte Carlo integration become vastly inef¬ 
ficient with increasing dimensionality. 

Therefore, we use importance sampling, a more gen¬ 
eral form of Monte Carlo integration, to calculate the in¬ 
tegral in Equation [2] an d make this problem c omputation- 
ally tractable. Following iFord &: Gregor^ ((2003), we sample 
from a distribution g{9) with a known normalization. We 
multiply the numerator and denominator of the integrand 
in Equation [ 2 ] by g{6), 

pid\M) = f (4) 

J gW 

While the value of p{d\M) has not changed, Equation [4] can 
be estimated by drawing N samples from g{0), 


p{4m) 


1 

N 


6ir^g(6) 


p{4ei,M)p{ei\M) 

gifii) 


(5) 


The key aspect of having importance sampling work 
efficiently is to pick an appropriate giO). Assuming our pa¬ 
rameter space contains one dominant posterior mode, we 
choose a multivariate normal with mean vector /2 and co- 
variance matrix E for g{9). For each model considered, we 
will estimate p, and E from the coplanar MCMC runs de¬ 
scribed in S|4l Our parameterization for g{0) is P, K, esincu, 
ecoscu, and lj + M for each planet, the system’s orbital incli¬ 
nation isya, and one aju for each observatory. Since we are 

only interested in computing ratios of p{d\M), the priors in 
zero-point offsets in the calculation will cancel out. 

One good strategy with importance sampling is to pick 
a g{9) that is heavier in the tails than p(d\9, M)p{9\M). 
This makes it easier to sample from low probability parts 
of the posterior distribution and prevents any samples from 
resulting in extremely large weights. However, the chance 
of sampling from the posterior mode will decrease with in¬ 
creasing dimensionality, which could ultimately lead to an 

estimate of p{d\M) that is not efficient. 

One way around this is to sample from g{9) within some 
truncated subspace, T. This new distribution griP) is pro¬ 
portional to g{9) inside T and renormalized to be a proper 
probability density. Equation can be rewritten as 

/xp(4aI)«1 ^ p{d\9,,M)m\M) ^ 

gr{di) 


where / is a factor that specifies what fraction of 
p{d,\9, M)p{9i\M) lies within T. We can estimate / with an 
MCMG sample. By counting what fraction of our posterior 
samples fell within T, Jmcmc, we can rearrange Equation 

[6] to give us p{d\M). 


p{d\M) 


1 

N X fmcMC 


E 

0i~9r{S) 


pi49i,M)p{9,\M) 

gT0i) 


(7) 


IGuoI (|2012ll and IWeinberg et al.l (l2013ll provide more de¬ 
tailed prescriptions and investigations of this method. 

There are two competing effects when choosing the size 
of our subspace T. If T is large (i.e. occupies nearly all of the 
posterior distribution), then Jmcmc approaches 1 and we 
return to our basic importance sampling algorithm. If T oc¬ 
cupies a much smaller region, then we are more likely to sam¬ 
ple from near the posterior mode, but /mcmc approaches 0, 

making it difficult to accurately estimate p{d\M). We must 
carefully choose a T that will provide a robust estimate of 

p{d\M).\G u3 (l2012l l found that for a three- and four-planet 
system, truncating the posterior distribution from — 2 ct to 
+2 (t was roughly optimal, where a is the standard devia¬ 
tion of each respective model parameter. We tested several 
subspace sizes and found ±1.5o' worked well for our problem. 

To draw from this distribution, we create a vector z 
whose components are independent draws from a standard 
normal A/'(0,1) truncated at —1.5a and -1-1.5cr for each of our 
vector components. Each 9i is generated as p-\-Az, where A 
is the Cholesky decomposition of E. We generate 10^ sam¬ 
ples from g{9), which provided a robust estimate of p{d\M) 
for each of our five competing models described in El Fig- 
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Number of Importance Samples [ xio® ] 


Figure 2. The estimate of log(p(d|Af)) (see Equation [7|l as a function of number of importance samples for five different models 
described in J4] and Table [2] Each panel shows the convergence of one model and are ranked from least probable (bottom) to most 
probable (top): Adsout (solid, blue), Adain (dashed, cyan), (dash-dotted, green), Af 3 in-|-~ (dotted, red), and Af 4 (solid, black). The 
difference in the minimum and maximum y-values of each panel is kept constant to better compare the variability in log(p(d|A4)) for 
each model we considered. 


Table 2. xlff, log evidence [log(p((i| Al)] , and Bayes factors computed for comparison of a finite set of models. 


Model 

Description 

xlff 

log(p((ilA4)) 

Bayes factor 

Ms 

Mi 

Al3in-|-~ 

Main 

Alsout 

Ad 4 + injected fifth 
all four planets 

Adsin + ^^120 day sinusoid 

3 inner-most planets 

3 outer-most planets 

743.8±|i5® 
768.3±|'g 
779.7±|-^ 
877.4±|'| 
1045.7±f;^ 

-1255.159 

-1248.139 

-1251.671 

-1273.519 

-1345.970 

~ 10-3 

~ 30 
~ 10® 
~103i 


The Bayes factors in right-most column in Table [2] correspond to the ratio of adjacent values of p{d\Ai) in the column immediately to 

the left. 


ure [5] shows how well the importance sampling algorithm 
converged for each model. 


4 COPLANAR MODELS 

Before attempting to relax the coplanarity constraint, we 
wish to assess the evideirce for all the GJ 876 planets, e in 
particular. 

We apply RUNDMC to our cumulative set of 367 RV ob¬ 
servations using a coplanar model but allowing for a sys- 
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Table 3. Estimates of the osculating orbital elements for all the known GJ 876 planets from self-consistent, coplanar dynamical fits. 


Parameter 

Planet d 

Planet c 

Planet b 

Planet e 

p [d] 

1 .yo / QQQQ23 

30.0784±g;00;5 

61.087±0;gii 

124.50 ±1.33 

K [ms“l] 

6.01 ± 0.30 

87.93 ±0.36 

213.94±g;37 

3-31±gii 

U.ii/lViQ r, - .j 

m M. 

7.49±0:« 

266.2±|;| 

845.3±li:g 

16.58±i;ff 

a [AU] 

0.02183911±0;0“| 

0.119±g:Ofo 

0.135989±o®gggi? 

0.218585 ±0.000026 

0.3514 ±0.0025 

e 

0.2531 ±0.0031 

0.0368±g:gg« 

0.031±g:g3g 

esinuj 

—0 044+®-^^^ 

U.U^^dlQ Q5g 

0.2276 ± 0.0045 

0.0345±g:ggiO 

0.008±g;gi3 

ecosoj 

_n nQ7_(_0.050 

u.uy/diQ Q4g 

U.liUOdlQ QQ42 

U.UlOUdlQ QQ20 

-0.013±g:gii 

i [degrees] 


53.19± 

1.39 

1.43 


LO [degrees] 

-140.90±^|6g|2 

115.96±i:i 

110.80±g|f 

129.56±3®31^ 

M [degrees] 

301.87±i|6«2^ 

-220.68±1|| 

-285.38±|lf 

-175.89±i?04f 

Lj M [degrees] 

161.71±^;^f 

-104.74±g-63 

-174.62±g;tg 

—45.05±fQ’^gy 


Estimates are computed using 15.9, 50, and 84.1 percentiles. For other details, refer to Section[4| 


Table 4. Systematic offset and jitter estimates based on a coplanar orbital model. 


Dataset 

Offset [ms 3 ] 

Jitter [ms ^ ] 

HIRES Carnegie (pre-upgrade) 

50.48±g:3| 

2.38±g;|3 

HIRES Carnegie (post-upgrade) 

52.20±g;i? 

HIRES California (post-upgrade) 

8.18 ±0.94 

6.82±g;g3 

ELODIE 

-1903.81±|f| 

18.56±|:73 

CORALIE 

-1864.40±g:i:| 

20.36±3;30 

HARPS 

-1337.64±g|| 

1.63±g:l 


Estimates are computed using 15.9, 50, and 84.1 percentiles. For other details, refer to Section IMl 


tematic orbital inclination. In particular, we consider five 
different models: one with just the outer-most three planets 
(Adsout), one with just the inner-most three planets (Adain), 
one with the inner-most three planets plus a ~ 124-day sinu¬ 
soid to mimic either a fourth planet on a circular orbit or a 
naive model for a stellar activity signal (A43in-i-~), one with 
all four planets {Mi), and one with a putative fifth planet 
(Ms). 

We calculate five fully marginalized likelihoods (four 
Bayes factors) and summarize the results in Table O Our 
methodology for computing these probabilities are described 
in iSH 

The Bayes factor for Msin/Msout is decisively 

favoring Adsin. This was expected since Main accounts 
for the three most significant RV signals in GJ 876. The 
Bayes factor for Adsinn —/Main is Despite the in¬ 

creased parameterization of Ad 3 in+~, this model with four 
signals is decisively favored over Adsin. The Bayes factor for 
Ad4/Ad3in-i-~ is ~30. The ratio of the evidence for a fully 
interacting four planet model Ad 4 relative to an interact¬ 
ing three planet model with a decoupled ~124-day signal 
Ad 3 in+~ is modest. Regardless, a model with four signals is 
overwhelmingly preferred over a model with three signals. 

Strictly considering the results above, we would not be 
able to select Ad 4 over Ad 3 in+~. Nevertheless, we do not ex¬ 
pect this ~ 124-day ca ndidate to be due to (an alias of) stel¬ 
lar magnetic activity. I Rivera et al.l ll2010li found that after 
fitting for the inner-three planets to the Carnegie RVs alone, 
the residual periodogram strongly peaked at the planet e’s 


period. The periodogram did not contain a strong signal near 
the stellar rotation period (see 43.31) . so we do not suspect 
the 124-day signal is an alias. 

For Ads, we initialized RUNDMC simulations with a hy¬ 
pothetical planet / with the following properties: Pf = 
15.04 d, Kf = 3ms"\ e/ = 0.007, ujf = 78.3°, and 
Mf = 159.8°. This set of initial conditions w as inspired by 
an ad ditional candidate signal uncovered by I Jenkins et al.l 
(l2014ll using a global periodogram method and RVs ex¬ 
tracted f rom HARPS spectra using the HARPS-TERRA 
pipeline (lAnglada-Escude fc ButlCT 120121 '). One significant 
difference is that we use an RV sign al much closer to th e 
noise level and HARPS RVs based on ICorreia et al.l ll2010ll . 
However, this periodicity is suspect since it is close to an 
integer fraction of the two most prominent signals in the 
system, c {PR2) and b (Pj,/4). This makes our choice of 
Ms SL good test for model comparison, since the strength of 
this signal competes with the extra parameterization needed 
to model it. After burning-in, we obtain the following esti¬ 
mates for Ms: Pf = 15.07±o: 29 d, Kf = 0.42 Toiti ms“^, 
and e/ = 0.16±o:i2- Due to the low RV amplitude, uif and 
Mf took on any value from 0° to 360° almost uniformly. We 
randomly draw 10,000 posterior samples from our Markov 
chains for each of these cases and report our values of Xe// in 
Table [21 along with the results of our Bayesian model com¬ 
parison tests. These Xe// values incorporate a penalty term 
depend ent on radial velocity jitter (Equation 5: [Nelson et al.l 
ll2C)14ar) ). 

We test the coplanar models for short-term orbital sta- 
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Figure 3. Mutual inclination distribution for planet’s d and 
c. The dashed red histogram shows the initial 10,000 posterior 
samples from Set 1 before any stability tests are performed. The 
solid green histogram shows a sample of 1,000 systems from Set 3 
that were stable for at least 10^ years. 


bility using the hybrid integrator of MERCURY (|Chambersl 
llQOlii i. A subsample of 1,000 models is integrated for 10® 
years. If at any point a planet collides with another body or 
the semi-major axis of any planet changes by more than 50% 
of its original value (i.e. ^ 0-5)) 

then the simulation stops and is tagged as being unstable. 
This is the default setup for our dynamical integrations un¬ 
less stated otherwise. 

For Alsout, we find 93% of the posterior samples are 
stable for the duration of the integration. The mode of insta¬ 
bility for the remaining 7% was | Aa/a| of planet e exceeding 
0.5. This typically happened when Ce > 0.14. For both Msin 
and A43in-i-~, all of the posterior samples are stable for the 
same integration timespan. For M 4 , 99% of the posterior 
samples are stable. For Ms, only 18% of the posterior sam¬ 
ples were stable. 77% of models had |Aa/a| of e exceeding 
0.5, and at least 5% of models involve a planetary collision. 
Some simulations had both instabilities occur. 

In M 4 , both the inner and outer Laplace pair were 
locked in a 2:1 mean-motion resonance, i.e. at least one res¬ 
onant argument was librating (for details on resonant argu¬ 
ments, see 321, for nearly all of our posterior samples. The 
Laplace argument was also librating in nearly all of our sam¬ 
ples. 

Based on these previous works and the results of i)3.3l 
and an we conclude that M 4 is the best model to gen¬ 
erate these observations to date. With a four-planet model 
decisively chosen, we report estimates of orbital parameters 
and instrumental properties based on a coplanar model in 
Tables |3] and |4] respectively. 

Before moving onto our three-dimensional model, we 
consider the significance of the relativistic precession rate 
chrei for planet d, since it is estimated to have a moderate 
eccentricity (cd ~0.1) and is close to its host star. While the 
rapid resonant interactions amongst the outer-three planets 
dominate their orbital evolution, the inner-most planet is 
mostly dynamically decoupled and undergoes secular per¬ 
turbations. A value of Wrei comparable to the secular preces¬ 


sion rate d>sec tends to aggravate an instability (iFord et al.l 
I 2 OOOII , while an curei th at is much more rapid than LOaec could 
quench secular effects (I Adams fc LaughlirJl2006alf3) . 

We approximate the precession rate 


(MM 

( ad ^ 


^ degrees 

\MeJ 

V0.05 AU/ 

' \ldayj 

century 


(8) 


where the subscript d refers to planet d and each parameter 
is measured in the un its of their respective normalization 
ll Jordan fc Bakoj|2008ll . We evaluate Equation [S] for 1,000 
four-planet models described and find Aei = 3.45±o:o3 de¬ 
grees per century. I Baluev! (1201 if) found that the effects of 
correlated noise results in a smaller estimate of Cd, though 
this does not change curei significantly. The typical ujsec for 
the same models was 88 ± 2 degrees per century. 

Ultimately, we wish to obtain a large sample of stable 
models but being careful not to neglect any mechanisms that 
could potentially stabilize a significant fraction of our mod¬ 
els. Therefore, relativistic precession will not be included in 
our final long-term integrations. 


5 RESULTS FOR 3D MODEL + LONG-TERM 
ORBITAL STABILITY 

The planet-planet interactions in GJ 876 are so strong that 
the physical planet masses must be used in the modeling 
process. Next, we allow RUNDMC to consider the full range of 
parameters associated with a three-dimensional (3-d) orbit 
for a four-planet model. To provide a frame of reference, we 
set fid = 0°. For our set of initial conditions for RUNDMC, 
we start with our posterior samples from M 4 and perturb 
each inclination and ascending node value (excluding fid) by 
adding a number drawn randomly and uniformly between - 
1° and 4-1°. RUNDMC burned-in for about 20,000 generations 
until a long-term steady state was reached. By modeling 
these extra variables, our parameter space grows from 32 
dimensions (4 planets x 5 orbital elements {P, K, e, w, M} + 
isys+ 6 offsets -I- 5 jitters) to 38. 

Additionally, we can obtain even more precise estimates 
by constraining the 3-d orbits of the planets from a direct 
analysis of the RV data and by demanding orbital stability. 
In the end, we report 1,000 solutions where the planetary 
system is dynamically stable for at least 10^ years. Our full 
procedure is described below and condensed into Table (5] 
We start by randomly drawing 10,000 posterior sam¬ 
ples from our initial 3-D MCMC run (Set 1) for our first set 
of stability tests. The model parameter estimates for these 
samples are shown in Tables [6] and 0 These systems are in¬ 
tegrated for 10® years and our conditions for stability are 
identical to those mentioned in 311 Only about 11% of our 
initial conditions passed our stability criterion; 48% of our 
sample show planet d falling into the central star, and 41% of 
our sample show planet e’s semi-major axis suddenly chang¬ 
ing by |Aa/a| > 0.5. We visually inspected the results and 
determined which model parameters were most important 
for distinguishing between “stable” and “unstable” regimes 
of parameter space. We find that the mutual inclination be¬ 
tween planet pairs is the most sensitive parameter in regards 
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Table 5. The different orbital parameter constraints applied to RUNDMC. 


Sample set 


Constraints during RUNDMC 


Then Integrate For... 

Set 1 

0° < >kdc < 180° 

0° < 'kcb < 180° 

0° < <I>be < 180° 

10®yr 

Set 2 

0° < <I>dc < 60° 

120° < <I>dc < 180° 

0° < 4>cb < 180° 

0° < <I>be < 13° 

10®yr 

Set 3 

0° < <I>dc < 60° 

120° < 4-dc < 180° 

0° < 4>cb < 4° 

0° < <I>be < 3^16 - $2,^" 

10’^ yr 


Here we show what range of parameter space RUNDMC is allowed to explore for each Set of simulations. The initial ensemble for a 
particular Set is generated based on stable regime found from the previous Set. 



0.88 0.90 0.92 0.94 0.96 0.98 1.00 

cos($(,J 


Figure 4. The joint cosine mutual inclination distribution for planets c and b (vertical axis) and b and e (horizontal axis). Contours 
map the approximate l-cr (68%) and 2-fT (95%) credible regions for the initial 10,000 posterior samples from Set 1 before any stability 
tests are performed. Horizontal hatch marks indicate the region for obtaining Set 2 and the vertical hatch marks indicate the region for 
obtaining Set 3, both described in f[5]and Table [5] Green dots show a sample of 1,000 systems from Set 3 that were stable for at least 
10^ years. 


to the system’s orbital stability. Henceforth, mutual inclina¬ 
tions between adjacent pairs of planets will be labeled as 
such: e.g. <f>dc for the mutual inclination between planets d 
and c (see Equation [T|). For planet d, 60° < i&dc < 120° 
causes d to undergo Kozai-like perturbations, pumping its 
eccentricity enough so that its periastron crosses the stel¬ 
lar surface (Figure [3|. Figure [4] shows the joint parameter 
distribution for cos(<l?cb) and cos('Fbe) based on Set 1 (con¬ 
tours). We choose this parameterization over <l?cb and i&be to 
visually demonstrate how consistent our solutions are with 
a coplanar system. For planet e, all initial conditions with 


i&be > 13° (cos(t&be) < 0.9781) led to a sudden and large 
change in planet e’s semi-major axis. Therefore, we perform 
a new set of RUN DMC simulations in which we restrict the pa¬ 
rameter space to exclude the unstable mutual inclinations 
found above, i.e. 60° < "Edc < 120°, "Fbe > 13°. Although 
we see some large values of 4>cb correspond to an instabil¬ 
ity, we do not restrict the range of <&cb for analysis in Set 2. 
RUNDMC is applied again on the same RV dataset, but now 
the resulting posterior sample (Set 2) is not as heavily di¬ 
luted by wildly unstable models. 

We repeat the above procedure one more time, inte- 
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Table 6. Estimates of the osculating orbital elements for all the known GJ 876 planets from self-consistent, 3-d dynamical fits. The top 
value in each cell are estimated directly from the radial velocities, and the bottom value imposes dynamical stability for 10^ years. 


Parameter 

Planet d 

Planet c 

Planet b 

Planet e 

p [d] 

1 QQ 7 QQ 1 _|_0.000037 
i.yc5royi±Q QQQ034 

1 Q‘^7S7n-{-0000025 
i.yo<o<uinQ QQQQ2g 

30.0758±0:g0f| 

30.0766±8:g073 

61.094 ±0.013 
61.087±0;0ii 

123.50±8|| 

124.72±1;26 

K [ms“^] 

a 1 o_|_0.39 

O.iOdio 35 

« 11 _{_0.35 
'^•-‘--*-^0.40 

88.74±0:|6 

88.33±0;4J 

213.63±0|0 

213.71±0;« 

3-i8±g:lS 

3.44±0:« 

0.37Mq r, - I 
^ M* 

7-03±3:|5 

6.90±g:|; 

274.2±i0:i 

267.9±f;| 

850.5±16:6 

848.5±8^;f 

15.45±i:|9 

17.16±2:i| 

a [AU] 

U.UZlocjy^Oito. 00000026 

U.U^lOOyOUHlQ QQQQQQ 20 

0 1‘^'^oau-l-0.000026 
U.iOOyO'illlQ QQQQ23 

0 1000022 
U.iOOyOOltQ QQQQjg 

0 218609+0-000034 
U.ZiODUyillQ QQQQ 3 J^ 

0 218589+0-000026 
u.zioooymQ QQQQ 3 Q 

0 ‘34QC:i-|_0.0031 

U.O^yOlllQ QQ25 

o.3518±8:8824 

e 

n 11 q_(_0.053 

U.iiCiito.osO 

o.io 8 ±g:gto 

u.zoozicq QQ3j^ 

0 2539+^-®®^^ 

u.zouyniQ QQ 34 

n 03714-0.0019 

U.UO^IiHq QQ29 

0 0365+0-0019 

U.UUDOiHq QQ3Q 

0.046±g;8« 

0 QO-I _|_0.027 
U.UOidiQ Q26 

esino; 

-o.o38±8:gti 

-0.034±8;8|§ 

0 22 n 8 -l -0 00 S’^ 

U.Z^UOitQ Q072 

0 2259+^-®®^^ 

U.ZZOyHlQ QQ 5 Q 

0 0345+0-0021 
U.UU40ilio.o034 

0 oo4o_|_0.0022 
U.UU^UHIq QQ3g 

0 01 Q-l-0.027 
U.UiyUiQ Q 2 g 

0 004+0-014 

U.UU4:I1Iq QJ 3 

e cos u) 

-0.090±8:0t? 

-o.o88±8:8« 

— 0 1 2363-^*^^^^ 
^•-‘-^'^0=^0.0101 

—0 1160+^-®®^^ 
U.liOUdlQ QQg 5 

— 0 01 ‘17-1-0.0024 

U.UlO/dlQ QQ21 

— 0 01 ‘I 7 -I-O OO 22 

U.UIU/IHq QQ2]^ 

— 0 025+0-033 

U.UZOilio.o40 

— 0 006+0-020 
U.UUDHIq Q3g 

i [degrees] 

87.44±«|7 

88 26+^®-^® 
OO.ZDIII 4 Q gQ 

51.63±3gl 

53.06±1;^ 

52.61±Hi 

52.82±l';|| 

55.51±5;gg 

53.29±8;8f 

oj [degrees] 

-139.25±||®9?° 

162.52±6;f7 

119.26±§:i 

117.12±i;™ 

111.81±|:fJ 

112.27±|;fl 

122.34±3g^005 

-54.2±23;8 

M [degrees] 

301.24±|^8*g2 

300.58±ii«|2 

-223.46±3:^ 

-221.76±f;88 

-286.01±|f8 

-286.82±|;2| 

-184.28±ilJ9 

—92 07-1-154.83 

u) -\- M [degrees] 

162.52±6'67 

162.28±6;^ 

-104.33±8'^f 

-104.60±8;|[ 

-174.26±8:S7 

-174.64±0:« 

-54.24±i;6| 

-42.46±ll9f2 

AQ [degrees] 


— 2 72+2-02 

AO L — '^^1.73 

9Q 1 0.96 
±.zyiiio.g9 


_c; 01 4-15.17 
AOi. — ^•^-‘-^10.91 

1 9q_l4.00 

i.zyit4 58 

4> [degrees] 


^ 6.20 

^ 2.60 


^ 28.5 

be < .j, gy 


Estimates are computed using 15.9, 50, and 84.1 percentiles. Upper limits on 4> are based on a 99% credible interval. For other details, 

refer to Section [5] 


Table 7. Systematic offset and jitter estimates based on a 3-d orbital model. 


Dataset 

Offset [ms 3 ] 

Jitter [ms ^ ] 

HIRES Carnegie (pre-upgrade) 

50.44 ±0.34 
50.49 ±0.31 

2.40±g:25 

HIRES Carnegie (post-upgrade) 

52.04±8;- 

52.12±8:65 

2.36±g;23 

HIRES California (post-upgrade) 

8.02 ± 0.90 
8.05±8:9i 

6.76±g:^l 

6.93±g:^t 

ELODIE 

-1904.25±||f 

-1903.83±|;|| 

18.97±|-go 

18.69±5;1? 

CORALIE 

-1864.32±8|8 

-1864.48±8;|| 

19.96±i:0i 

19-94±iii 

HARPS 

-1337.68±gi8 

-1337.58±g;|i 

1 K2+0-27 
i.JZHiQ 24 

1 c7_(_0.25 
-*-•^'^0.22 


Estimates are computed using 15.9, 50, and 84.1 percentiles. For other details, refer to Section 13.21 
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grating the Set 2 for 10® years. We find the initial condi¬ 
tions are unstable, unless 0° < "hcb < 4° and 0° < "hbe < 
3^16 — . We run RUNDMC once more while implement¬ 

ing these tighter constraints to obtain a new set of posterior 
samples for Set 3, which are subsequently integrated for 10^ 
years. Parameter estimates of 1,000 stable solutions from 
Set 3 are shown as the bottom value of each cell in Table [6l 
as the solid green histogram in Figure [3] and as green dots 
in Figured! 

Performing everything described above, we find the ra¬ 
dial velocities alone place an upper limit of "Feb < 6.20° and 
$be < 28.5° based on a 99% credible interval. However if we 
expect the system to remain stable over the course of 10^ 
years, the system must be roughly coplanar: tf>cb < 2.60° 
and 4?be < 7.87° based on a 99% credible interval. 


6 BEHAVIOR OF RESONANT ANGLES 
6.1 Eccentricity Resonances 

With the final 3-d orbital solutions from Set 3 of US) we in¬ 
vestigate the behavior of the critical angles associated with 
the mean-motion resonances relevant for this system. For 
both our stable coplanar and 3-d orbital models, we com¬ 
pute the root-mean-square of the variability in each angle 
xy^. For a system undergoing small amplitude sinusoidal 
libration, this is an excellent approximation for the libration 
amplitude. 

For the c and b pair, the angles associated with the 2:1 
MMR are the resonant angles, = Ac — 2Xb -I- zuc and 
(j)f‘ = Ac — 2At -I- zjb, and the secular angle zuc — zuh, where 
A and zu are the mean longitude and longitude of pericenter 
respectively. We find each of these angles is librating with 
low amplitude about 0°. Figure [5] shows the distribution of 
the libration amplitude for angles associated with the c and b 
pair. Similarly for the b and e pair, the angles associated with 
the 2:1 MMR are the resonant angles, Rlp = At — 2Ae -I- rot 
and (^c° = At — 2Ae -l-TOe, and the secular angle rot — zu^- We 
find librates about about 0° and the other two angles are 
circulating. Figure [5] shows the distribution of the libration 
amplitude for angles associated with the b and e pair. For 
the c and e pair, the angles associated with the 4:1 MMR 
are the resonant angles, = Ac — 4Ae -|- 3zuc, (pT = Ac — 

4Ae ‘2zUc 4- ZUe^ (p2 ~ Ac — 4Ae 4- ZUc 4- 2zUe^ and 03 = 

Ac — 4 Ae-|- 3 n 7 e, and the secular angle zuc — zus- To distinguish 
these angles, the subscript refers to the multiplier in front 
of zus- We find that four of the five angles are circulating. 
00 ° librates about 0° with low to medium amplitude. Figure 
[7] shows the distribution of the libration amplitude for 0o°. 
We report all of our libration amplitudes for coplanar and 
dynamically stable 3-d orbital models in Table [ 8 | 

We find the secular angle oJb — uJe is circulating, in 
contrast to previous studies that reported libration about 
180°. This underscores the importance of performing self- 
consistent dynamical and statistical analyses when charac¬ 
terizing the evolution of interacting planetary systems. 

The measured amplitude is interesting as it is neither 
so small as to imply strong dissipation nor so large as to 
suggest the absence of damping. Given the importance of 
this result, we performed additional tests to verify that our 


Table 8. Libration Amplitudes of Resonant Angles. 


Angle 

Amplitude, coplanar [°] 

Amplitude, 3D [°] 

pf 

4.0±1-6 

7-o±l:! 

Pt 

13.4±3-2 

18.69±i:5 

TJJ Q 

14.7±|:g 

20.6±§:| 

Pf 

PT 

31.3±ii7® 

47.9±y 

circulating 

circulating 

ZUi) — ZUe. 

circulating 

circulating 

00° 

67.5±i:| 

103.24i!|J 

0f 

circulating 

circulating 

0f 

circulating 

circulating 

PT 

circulating 

circulating 


circulating 

circulating 

4^Laplace 

33.0±1234 

50.5±l5% 


See Section!^ for details. 



Figure 6. Libration amplitude distributions for only librating 
angle of the b and e planet pair in a 2:1 mean-motion resonance. 
We compare libration amplitudes for a coplanar (dashed blue) 
and 3-d (solid orange) orbital model. For the vast majority of our 
sets of initial conditions, the other resonant argument and the 
secular angle are circulating. 


algorithm accurately characterizes the libration amplitude 
of the Laplace angle for both regular and chaotic systems. 
The methods and more detailed results of these test are 
presented in Appendix El 

Upon the d iscovery of the outer-most planet, 
I Rivera et al. found that their best fit systems exhibit 

a 4:2:1 Laplace resonance. The associated angle, pLapiacs = 
Ac — 3A6 4- 2Ae, evolves chaotically. Based on the Carnegie 
RVs alone, they found that pLapiace librates about 0° with 
an amplitude of 40 ± 13° for a coplanar four-planet model. 
Figure [5] shows our results for the posterior distribution for 
the libration amplitude for the Laplace argument for both 
the coplanar (33±9^3^°) and 3-d (50.74ri5®°) cases. The mea¬ 
sured amplitude is interesting as it is neither so small as to 
imply strong dissipation nor so large as to suggest the ab¬ 
sence of damping. Given the importance of this result, we 
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Figure 5. Libration amplitude distributions for the three resonant angles of the c and b planet pair in a 2:1 mean-motion resonance. 
We compare libration amplitudes for a coplanar (dashed blue) and 3-d (solid orange) orbital model. The two angles associated with the 
2:1 mean motion resonance (top and middle panels) are librating about 0° with low amplitude. The secular angle (bottom panel) also 
librates about 0° with low amplitude. 


performed additional tests to verify that our algorithm accu¬ 
rately characterizes the libration amplitude of the Laplace 
angle for both regular and chaotic systems. The methods 
and more detailed results of these test are presented in Ap¬ 
pendix 

6.2 Inclination Resonances 

Several studies have looked at inclination excitation and 


inclination resonance capture during planet migration 

phases (iThommes & Lissaueil l2003l: iLee & Thommes 

l2009l; 

iLibert & TsieanisI 20091: Tevssandier & Terauemll2014 

). The 


general conclusion in these studies is that inclination reso¬ 
nances can form during the form ation phase bu t are li kely to 
be short-lived. A recent study bv iBarnes et al.l (l2015h inves¬ 
tigated the orbital evolution of MMR systems with mutually 
inclined orbits in the post-formation phase. They find that 
their synthetic systems evolve chaotically but still retain dy¬ 
namical stability and a MMR. Drawing upon these conclu¬ 
sions, they model several real RV systems near a MMR to 
look for similar behavior. They find that HD 73526 (a 2:1 


MMR system) and HD 60532 (a 3:1 MMR system) could be 
evolving chaotically given a set of initial conditions listed in 
their Tables 9 and 10 respectively. However, the inclination 
resonance arguments associated with these systems were all 
circulating, indicating that they were not in an inclination 
resonance. To date, no exoplanet system has a detected in¬ 
clination resonance. 

We check for an inclination resonance in GJ 876 with 
our dynamically stable posterior samples. First, we trans¬ 
form our set of initial conditions into the invariable plancj3 
then integrate these for 1,000 years, logging the output ev¬ 
ery 10 days. We compute the angles associated with the 
inclination-resonance for the inner pair 4>inci = 4A{, — 2Ac — 
Dc — fli, and the outer pair ~ 4Ae — 2At — Dt — De. Visu¬ 

ally inspecting the orbital evolution over 100 years of these 
models, we find ipind s-nd (pi^d both circulate on timescales 
of ~5 years. While it is possible that pind may ocassionally 
librate for a short period of time, this is sufficiently rare that 


^ https://github.com/RoryBarnes/InvPlane 
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Figure 7. Libration amplitude distribution for the only librating 
angle of the c and e planet pair relative to a 4:1 mean-motion res¬ 
onance. We compare libration amplitudes for a coplanar (dashed 
blue) and 3-d (solid orange) orbital model. For the vast majority 
of our sets of initial conditions, the three other resonant argu¬ 
ments and the secular angle are circulating. 



Figure 8. Libration amplitude distributions for the Laplace ar¬ 
gument of the c, b, and e resonant trio. We compare libration 
amplitudes for a coplanar (dashed blue) and 3-d (solid orange) 
orbital model. 


we did not observe any such events over this baseline for all 
of our models. 


7 SUMMARY AND DISCUSSION 

We have meaningfully constrained the three-dimensional or¬ 
bital architecture of the GJ 876 planetary system, based on 
367 Doppler observations from several observing sites. 

To verify the nature of the ~ 124-day signal, we per¬ 
formed Bayesian model comparison of five different physical 
models, spanning three to five planets. Since each evalua¬ 
tion of the likelihood requires an n-body integration for a 
strongly interacting planetary system like GJ 876, Bayesian 


model comparison becomes computationally costly. There¬ 
fore, it was particularly important that we apply efficient 
and parallelizable algorithm. We refined and applied a 
modified importance sampling algorithm to compute the 
fully marginalized likelihood, or Bayesian evidence, start¬ 
ing from a posterior sample computed via MGMC meth - 
ods jFord fc Gregorvll2007l : lGuol2012l : IWeinberg et al.|[201^ ). 

The algorithm parallelizes readily and was implemented o n 
GPUs using the Swarm-NG framework llDindar et akluMS ). 
While previous studies have computed B ayes factors using 
other algorithms (e.g. nested sam pling dFeroz et al.l I 2 OO 9 I : 
iKippin j gOlj^ JPlacelc et ^ [20l3) , geometric-path Monte 
Garlo ( Hou et al.l 2014j) ). these studies have assumed that 
the motion can be described as the linear superposition of 
Keplerian orbits, which is unsuitable for strongly interact¬ 
ing planetary systems such as GJ 876. We believe this study 
to be the first example of rigorous Bayesian model compar¬ 
ison applied to strongly interacting planetary systems. This 
algorithm is relatively easy to implement and worked well 
even for our high-dimensional (~30-40 parameter) models. 


We determined that a four-planet model is most ap¬ 
propriate for the present data, based on a self-consistent 
Bayesian and n-body analysis ( i)3.3l and 13.411 . When com¬ 
puting the evidence for our finite set of models, we can de¬ 
cisively choose a model with four signals with Bayes factors 
exceeding 10®. We find a Bayes factor of ~30 when compar¬ 
ing a four planet model {Mi) to a model with three planets 
plus a decoupled sinusoid, which could result from either a 
fourth planet in a circular orbit or a stellar activity signal 

(A43in+~). 

On one hand, a Bayes factor of ^30 is not large enough 
by itself to definitively choose the four planet model. On the 
other hand, we find no reason to suspect that stellar activity 
is masquerading as a planetary signal. Looking at activity- 
sensitive indicators in publicly available HARPS spectra, 
we measure a stellar rotation period of 95 ± 1 days. This 
does not directly correspond with any of the planets’ orbital 
periods. As expected, the stellar rotation period and annual 
observing cycle lead to an alias with a frequency of ~128 
days. In our coplanar model, the orbital period of e is 124.5± 
1.3, differing from the 128-day alias at slightly more than 
2 -(t level. If the signal were due solely to aliasing of the 
annual observing cycle with the stellar rotation period, then 
there should be an even stronger signal corresponding to the 
stellar rotation in the periodogram of the RV residuals (after 
subtracting the RV sig nal due to th e inner-three planets). 
I Rivera et al.l ll2010l ) and iBaluevI ll201ll ) looked for but did not 
find this effect, concluding that the signal is best explained 
by planet e. 

Through numerical integration s of systems where p lanet 
e was treated as a test particle, iRivera et al.l ll20ld) ex¬ 
plored the parameter space near the best-fit GJ 876 solu¬ 
tion to the Garnegie RVs. These tests suggested that the 
true system might be chaotic and was likely surrounded by 
regions of phase s pace with short-lived (unstable) orbits. 
iMartf et al.l (l2013l) studied the two reported 4-planet solu¬ 
tions to_thej;adial_velocity_data using the MEGNO method 
(e.g. iGincotta fc Simdl (120001) ) and found that these trajec¬ 
tories were both chaotic and in the Laplace resonance. We 
show that the most probable Lyapunov time to be ~10 years 
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for both our coplanar and dynamically stable 3-d models 
f AppendixIASII. This tim escale is consistent with the results 
of iBatygin et alj (l2015fl which were based on combining a 
simplihed dynamical model and a point estimate for orbital 
parameters. 

Strongly interacting planetary systems like GJ 876 
could have unusually large transit timing variation am¬ 
plitudes, enabling detailed characterization via the transit 
timing technique. Upcoming missions such as TESS and 
PLATO would easily detect and confirm a similar system 
of transiting planets, but a detailed interpretation could be 
challenging given the anticipated observing timespans. In 
particular, the perturbations f rom GJ 876 e become e vident 
only on multi-year timescales dLibert fc Renneill2013f l. 

While previous studies have assumed that the planets 
in GJ 876 follow coplanar orbits, we relax the assumption 
of coplanarity and allow for non-coplanar orbital configura¬ 
tions. We analyze the Doppler observations of GJ 876 and 
hnd that the three planets participating in the Laplace reso¬ 
nance must be nearly coplanar. The 99% credible upper lim¬ 
its on the mutual inclination constrained by just the RVs is 
impressive for both planet pairs: "Lcb < 6.20° for the c and 
h pair and c&be < 28.5° for the b and e pair. Demanding 
orbital stability further restricts this range providing pre¬ 
cise constraints on the mutual inclinations: "hcb < 2.60° and 
"Lbe < 7.87°. A plot of the posterior samples suggest that the 
orbits of these planets are consistent with a coplanar system 
(Figure SI. Despite its rather unique orbital architecture, it 
seems that GJ 876 fits in with a pop ulation of M dwarf 
syste ms with several, coplanar planets (iBallard fc Johnsorj 
|2014) . 

By performing the hrst self-consistent, Bayesian analy¬ 
sis of the four planets in GJ 876, we are able to accurately 
characterize the current dynamical state of all four planets 
and particularly the evolution of the resonant angle associ¬ 
ated with the Laplace resonance. We measure the amplitude 
of variations to be degrees for coplanar models or 

50.5±Io% degrees for fully three-dimensional models. 

When measuring a positive definite quantity, obser¬ 
vational uncertainties can bias measurements, particularly 
when using point estimates. For example, the best-fit mod¬ 
els of Doppler observations of a population of planets with 
nearly circular orbits will typically overestimate the plan¬ 
ets’ orbital eccentricities JZakamska et al.l 1201 l| j. particu¬ 
larly when the Doppler amplitude is only a factor of a few 
greater than the measurement precision. This motivated our 
Bayesian approach to characterizing the amplitude of varia¬ 
tions of the Laplace angle. Additionally, we performed tests 
of our algorithms using simulated planetary systems that 
confirm our algorithm accurately characterizes the behavior 
of the Laplace angle f Appendix I Alll . We conclude that for¬ 
mation theories for the GJ 876 system need to explain not 
only the resonant structure, but also the chaotic evolution 
of the Laplace angle and its sizable amplitude. 

The near integer ratio in the planets’ orbital periods 
are not likely a result of happenstance. The probability of 
in situ formation yielding such a system is further reduced 
by the need to become trapped in a chaotic, but long-lived, 
multi-body resonance. Within the context of current planet 
formation models, the system most likely reached its cur¬ 


rent resonant configuration as the result of disk migration. 
Migration through a smooth gas disk would be expected to 
result in strong damping of eccentricities, driving the sys¬ 
tem to a state where the resonant angle librates regularly 
with small amplitude. Many of our simulations of a smooth 
migration with eccentricity damping lead to an amplitude of 
~1°, much smaller than we measure for the current system 
configuration. 

While the observed large libration amplitude and 
chaotic evolution of resonant angles is c ontrary to the predic - 
tions of the simplest migration models. [Batygin et al.l (l2015lj 
recently proposed that that the observed libration of the 
resonant angles can be explained purely by a turbulent mi¬ 
gration. We propose an alternative formation mechanism 
based on a phase of smooth disk migration which terminates 
abruptly (Appendix|^. In our simplistic migration models, 
turning off the eccentricity damping impulsively caused the 
libration amplitude for the Laplace angle to rise from ~ 1° 
to ~tens of degrees. Such a scenario could arise naturally as 
a result of rapid dispersal of inner disk via photoevaporation. 

Another potential possible formation scenario involves 
migration through a gas disk trapping the planets in reso¬ 
nance, followed by a phase of p lanet or planetessimal scat¬ 
tering. ([Chatteriee fc Fordll2015l l showed that planetessimal 
scattering will naturally drive a planetary system initially 
in a 2:1 mean motion resonance farther apart, exciting ec¬ 
centricities and perhaps even breaking the resonance. 

It is p ossible for MMRs to ari se from pure planet-planet 
scattering. [Raymond et al.l ll2008lf found that MMRs could 
form through the ejection of one planet, based on simula¬ 
tions starting with three planets with drastically different 
masses. One can distinguish between MMRs formed through 
scattering and migration mechanisms by the mass and or¬ 
bital properties of the observed planets. Scattering usually 
yields planet pairs in higher order resonances (second or 
greater), with larger semi-major axes (>1AU), and/or with 
a more massive outer companion. The latter result could 
explain the c and b pair, but after considering e and the ob¬ 
served Laplace resonance, planet-planet scattering alone be¬ 
comes a less likely formation channe l for the GJ 876 planets. 
Indeed. [Moeckel fc Armitaed (|2012|j run a set of simulations 
modeling multi-planet systems during disk clearing phase 
and subsequent gas-free (purely Newtonian) phase and find 
that systems were much more likely to form resonant chains 
if scattering events did not occur. In the presence of a plan¬ 
etessimal disk, alread y “marginally stable” s ystems can also 
form resonant chains jRavmond et al.ll2009l j. 

We encourage future studies to explore the predictions 
of these formation models for comparison to the GJ 876 
system. Continued long-term RV monitoring and/or astro- 
metric observations from GAIA could continue to improve 
the dynamical constraints on this landmark system. 


ACKNOWLEDGEMENTS 

We would like to thank our referee, Edward Thommes, for 
his constructive feedback that improved the manuscript. We 
would like to thank Geoff Marcy and the entire of the Cali¬ 
fornia Planet Survey team for their long-term commitment 


© 0000 RAS, MNRAS 000, 000-000 






























3D Resonance in GJ 876 15 


to high-precision RVs for the GJ 876 system. B.E.N. would 
like to thank Phil Gregory and Tom Loredo for useful sug¬ 
gestions regarding our methodology and presentation of our 
importance sampling algorithm. Additionally, he would also 
like to thank Rory Barnes and Russell Deitrick for conver¬ 
sations regarding chaos and inclination resonances in multi¬ 
planet systems and Roman Baluev for useful feedback re¬ 
garding the significance of correlated noise. 

P.R. acknowledges support from NSF grant AST- 
1126413 and the Center for Exoplanets and Habitable 
Worlds. M.J.P. gratefully acknowledges the NASA Origins 
of Solar Systems Program grant NNX13A124G. E.B.F. and 
J.T.W. acknowledge NASA Keck PI Data Awards, admin¬ 
istered by the NASA Exoplanet Science Institute, including 
awards 2007B N095Hr, 2010A N147Hr, 2011A&B N141Hr, 
& 2012A N129Hr. This research was supported by NASA 
Origins of Solar Systems grant NNX09AB35G. The authors 
acknowledge the University of Florida High Performance 
Computing Center and the Pennsylvania State Research 
Computing and Advanced Cyberinfrastructure Group for 
providing computational resources and support that have 
contributed to the results reported within this paper. The 
Center for Exoplanets and Habitable Worlds is supported 
by the Pennsylvania State University, the Eberly College 
of Science, and the Pennsylvania Space Grant Consortium. 
We extend special thanks to those of Hawai’ian ancestry 
on whose sacred mountain of Manna Kea we are privileged 
to be guests. Without their generous hospitality, the Keck 
observations presented herein would not have been possible. 


REFERENCES 

Adams, F. C. & Laughlin, G. 2006a, ApJ, 649, 992 
Adams, F. C. & Laughlin, G. 2006b, International Journal 
of Modern Physics D, 15, 2133 
Anglada-Escude, G. & Butler, R. P. 2012, ApJS, 200, 15 
Anglada-Escude, G., Lopez-Morales, M., & Chambers, 
J. E. 2010, ApJ, 709, 168 
Ballard, S. & Johnson, J. A. 2014, ArXiv e-prints 
Baluev, R. V. 2011, Celestial Mechanics and Dynamical 
Astronomy, 111, 235 

Barnes, R., Deitrick, R., Greenberg, R., Quinn, T. R., & 
Raymond, S. N. 2015, ApJ, 801, 101 
Batygin, K., Deck, K. M., & Holman, M. J. 2015, AJ, 149, 
167 

Batygin, K. & Morbidelli, A. 2013, AJ, 145, 1 
Bean, J. L. & Seifahrt, A. 2009, A&A, 496, 249 
Beauge, C., Ferraz-Mello, S., & Michtchenko, T. A. 2003, 
ApJ, 593, 1124 

Benedict, G. F., McArthur, B. E., Forveille, T., et al. 2002, 
ApJ, 581, L115 

Boisse, I., Bouchy, F., Hebrard, G., et al. 2011, A&A, 528, 
A4 

Butler, R. P., Marcy, G. W., Williams, E., et al. 1996, 
PASP, 108, 500 

Chambers, J. E. 1999, MNRAS, 304, 793 
Chatterjee, S. & Ford, E. B. 2015, ApJ, 803, 33 
Cincotta, P. M. & Simo, C. 2000, A&AS, 147, 205 


Correia, A. C. M., Couetdic, J., Laskar, J., et al. 2010, 
A&A, 511, A21 

Delfosse, X., Forveille, T., Mayor, M., et al. 1998, A&A, 
338, L67 

Dindar, S., Ford, E. B., Juric, M., et al. 2013, New Astron., 
23, 6 

Dumusque, X., Boisse, I., & Santos, N. C. 2014, ApJ, 796, 
132 

Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 
398, 1601 

Ford, E. B. & Gregory, P. C. 2007, in Astronomical Soci¬ 
ety of the Pacific Conference Series, Vol. 371, Statistical 
Challenges in Modern Astronomy IV, ed. G. J. Babu & 
E. D. Feigelson, 189—h 

Ford, E. B., Joshi, K. J., Rasio, F. A., & Zbarsky, B. 2000, 
ApJ, 528, 336 

Ford, E. B., Moorhead, A. V., & Veras, D. 2011, 
Bayesian Anal., 6, 475 

France, K., Froning, C. S., Linsky, J. L., et al. 2013, ApJ, 
763, 149 

France, K., Linsky, J. L., Tian, F., Froning, C. S., & 
Roberge, A. 2012, ApJ, 750, L32 
Gerlach, E. & Haghighipour, N. 2012, Celestial Mechanics 
and Dynamical Astronomy, 113, 35 
Gozdziewski, K., Bois, E., & Maciejewski, A. J. 2002, MN¬ 
RAS, 332, 839 

Guo, P.-C. 2012, PhD thesis. University of Florida 
Haghighipour, N., Couetdic, J., Varadi, F., & Moore, W. B. 
2003, ApJ, 596, 1332 

Hon, F., Goodman, J., & Hogg, D. W. 2014, ArXiv e-prints 
Jenkins, J. S., Yoma, N. B., Rojo, P., Mahu, R., & Wuth, 
J. 2014, MNRAS, 441, 2253 
Ji, J., Li, G., & Liu, L. 2002, ApJ, 572, 1041 
Ji, J. H. & Liu, L. 2006, Acta Astronomica Sinica, 47, 402 
Johnson, J. A., Payne, M., Howard, A. W., et al. 2011, AJ, 
141, 16 

Jones, B. W., Sleep, P. N., & Chambers, J. E. 2001, A&A, 
366, 254 

Jordan, A. & Bakos, G. A. 2008, ApJ, 685, 543 
Kammer, J. A., Knutson, H. A., Howard, A. W., et al. 2014, 
ApJ, 781, 103 

Kinoshita, H. & Nakai, H. 2001, Publ. Astron. Soc. Japan, 
53, L25 

Kipping, D. M. 2013, MNRAS, 434, L51 
Kley, W., Lee, M. H., Murray, N., & Peale, S. J. 2005, 
A&A, 437, 727 

Kley, W., Peitz, J., & Bryden, G. 2004, A&A, 414, 735 
Kokubo, E., Yoshinaga, K., & Makino, J. 1998, MNRAS, 
297, 1067 

Laughlin, G., Butler, R. P., Fischer, D. A., et al. 2005, ApJ, 
622, 1182 

Laughlin, G. & Chambers, J. E. 2001, ApJ, 551, L109 
Lee, M. H. 2004, ApJ, 611, 517 
Lee, M. H. & Peale, S. J. 2002, ApJ, 567, 596 
Lee, M. H. & Thommes, E. W. 2009, ApJ, 702, 1662 
Lega, E., Morbidelli, A., & Nesvorny, D. 2013, MNRAS, 
431, 3494 

Leinert, C., Henry, T., Glindemann, A., & McCarthy, Jr., 
D. W. 1997, A&A, 325, 159 

Libert, A.-S. & Renner, S. 2013, MNRAS, 430, 1369 


© 0000 RAS, MNRAS 000, 000-000 


16 B. Nelson et al. 


Libert, A.-S. & Tsiganis, K. 2009, A&A, 493, 677 
Lichtenberg, A. J. & Lieberman, M. A. 1992, Regular and 
chaotic dynamics (Springer) 

Luhman, K. L. & Jayawardhana, R. 2002, ApJ, 566, 1132 
Marcy, G. W. & Butler, R. P. 1992, PASP, 104, 270 
Marcy, G. W., Butler, R. P., Fischer, D., et al. 2001, ApJ, 
556, 296 

Marcy, G. W., Butler, R. P., Vogt, S. S., Fischer, D., & 
Lissauer, J. J. 1998, ApJ, 505, L147 
Marti, J. G., Giuppone, G. A., & Beauge, C. 2013, MNRAS, 
433, 928 

Milani, A. & Nobili, A. M. 1983, Gelestial Mechanics, 31, 
213 

Moeckel, N. & Armitage, P. J. 2012, MNRAS, 419, 366 
Montet, B. T., Grepp, J. R., Johnson, J. A., Howard, A. W., 
& Marcy, G. W. 2014, ApJ, 781, 28 
Murray, N., Paskowitz, M., & Holman, M. 2002, ApJ, 565, 
608 

Nelson, B. E., Ford, E. B., & Payne, M. J. 2014a, ApJS, 

210, 11 

Nelson, B. E., Ford, E. B., Wright, J. T., et al. 2014b, 
MNRAS, 441, 442 

Patience, J., White, R. J., Ghez, A. M., et al. 2002, ApJ, 
581, 654 

Placek, B., Knuth, K. H., & Angerhausen, D. 2014, ApJ, 
795, 112 

Podlewska-Gaca, E., Papaloizou, J. C. B., & Szuszkiewicz, 
E. 2012, MNRAS, 421, 1736 
Rauch, K. P. & Holman, M. 1999, AJ, 117, 1087 
Raymond, S. N., Armitage, P. J., & Gorelick, N. 2009, ApJ, 
699, L88 

Raymond, S. N., Barnes, R., Armitage, P. J., & Gorelick, 
N. 2008, ApJ, 687, L107 

Rivera, E. & Haghighipour, N. 2007, MNRAS, 374, 599 
Rivera, E. J., Laughlin, G., Butler, R. P., et al. 2010, ApJ, 
719, 890 

Rivera, E. J. & Lissauer, J. J. 2001, ApJ, 558, 392 
Rivera, E. J., Lissauer, J. J., Butler, R. P., et al. 2005, ApJ, 
634, 625 

Robertson, P. & Mahadevan, S. 2014, ApJ, 793, L24 
Robertson, P., Mahadevan, S., Endl, M., & Roy, A. 2014, 
Science, 345, 440 

Shankland, P. D., Rivera, E. J., Laughlin, G., et al. 2006, 
ApJ, 653, 700 

Snellgrove, M. D., Papaloizou, J. C. B., & Nelson, R. P. 
2001, A&A, 374, 1092 

Tan, X., Payne, M. J., Lee, M. H., et al. 2013, ApJ, 777, 
101 

Teyssandier, J. & Terquem, C. 2014, MNRAS, 443, 568 
Thommes, E. W. 2005, ApJ, 626, 1033 
Thommes, E. W. & Lissauer, J. J. 2003, ApJ, 597, 566 
Thommes, E. W., Matsumura, S., & Rasio, F. A. 2008, 
Science, 321, 814 

Valenti, J. A., Butler, R. P., & Marcy, G. W. 1995, PASP, 
107, 966 

van Leeuwen, F., ed. 2007, Astrophysics and Space Science 
Library, Vol. 350, Hipparcos, the New Reduction of the 
Raw Data 

Veras, D. 2007, Celestial Mechanics and Dynamical Astron¬ 
omy, 99, 197 


Veras, D. & Ford, E. B. 2010, ApJ, 715, 803 
Vogt, S. S., Allen, S. L., Bigelow, B. C., et al. 1994, in 
Proc. SPIE Instrumentation in Astronomy VHI, David L. 
Grawford; Eric R. Craine; Eds., Vol. 2198, p. 362 
von Braun, K., Boyajian, T. S., van Belle, G. T., et al. 
2014, MNRAS, 438, 2413 

Weinberg, M. D., Yoon, L, & Katz, N. 2013, ArXiv e-prints 
Wisdom, J. 2006, AJ, 131, 2294 
Wisdom, J. & Holman, M. 1991, AJ, 102, 1528 
Wisdom, J. & Holman, M. 1992, AJ, 104, 2022 
Wisdom, J., Holman, M., & Touma, J. 1996, Fields Insti¬ 
tute Gommunications, Vol. 10, p. 217, 10, 217 
Zakamska, N. L., Pan, M., & Ford, E. B. 2011, MNRAS, 
410, 1895 

Zechmeister, M. & Kiirster, M. 2009, A&A, 496, 577 
Zhou, J.-L. & Sun, Y.-S. 2003, ApJ, 598, 1290 


© 0000 RAS, MNRAS 000, 000-000 


3D Resonance in GJ 876 17 


APPENDIX A: FITTING SYNTHETIC 
SYSTEMS WITH KNOWN PROPERTIES 

iBatvgin et alJ ll2015fl provide an elegant explanation of how 
such a chaotically evolving Laplace angle could have formed, 
and show that it can provide limits on the system’s early for¬ 
mation: smooth migration can only form systems with rel¬ 
atively small libration amplitude, and regular non-chaotic 
evolution; stochastic migration is required to form chaoti¬ 
cally evolving systems with large libration amplitudes. 

At face value, our results provide evidence for 
chaotic evolution of the Laplace argument in the GJ 876 sys¬ 
tem. Given the importance of this result, we performed ad¬ 
ditional tests to determine whether low-libration amplitude 
is likely to be erroneously classified chaotic orbital evolu¬ 
tion. We simulated RV observations of a variety of synthetic 
coplanar systems and then applied our RUN DMC algorithm to 
verify that the recovered parameter estimates (in particular 
the libration amplitude of the Laplace angle) are consistent 
with their input values. 


A1 Methodology: Creating Synthetic Planetary 
Systems and Datasets 

Batygin {private communication) provided orbital elements 
for the results of one of their stochastic migration simula¬ 
tions. We integrated this system forward, confirming that 
the evolution of the Laplace resonance angle, <j>Lapiace, is 
indeed chaotic, and has a libration amplitude ~80°. This 
synthetic system will be referred to as Sc^iaosiso henceforth 
(Figure [All left!. 

To complement this, we then performed a number of 
smooth migration simulations in which the resultant simu¬ 
lations are regular and have low libration amplitudes. These 
migration simulations used th e damped a and e method de¬ 
scribed in I Lee fc Peald (l2002l ), which was impl emented in a 
modified version of MERCURY llChambersIIl999l l. We set up 
3-planet systems so that each initial planetary period ratio 
was a little larger than 2:1 with their inner neighbor. Then, 
we applied the damped-migration model to the outer planet 
until it caught into a 2:1 resonance with the middle planet. 
We continued the forced, damped migration until this outer 
resonant pair caught into a 4:2:1 resonance with the inner¬ 
most planet. Finally, we removed all forms of damping from 
the system and let them evolve for ~10® years, ensuring that 
they were in an undamped equilibrium state. 

We performed a number of these simulations and se¬ 
lected a system from these with final (i.e. after turning off 
orbital damping) libration amplitude of ~20°. Most sim¬ 
ulations resulted in a Laplace argument with very small 
libration amplitude (~1°) during the damping phase and 
rose to several tens of degrees after the damping was re¬ 
moved. This particular system had a libration amplitude 
significantly smaller than the best-fit from allowing us 
to test whether our fitting procedure artificially increases the 
libration amplitude. This synthetic system will be referred 
to as S_Reg ;20 henceforth (Figureleft). 

We generate an RV dataset for both Schaos-.so and 
Si?eg :20 using the same RV time series with their associated 
offsets and jitters described in 1)3.21 We include the inner¬ 


most, non-resonant planet with the following orbital proper¬ 
ties and coplanar with the migrated planets: Pd = 1.937821 
days, md = 2.21 x 10“®Mq, Cd = 0.072, uJd = —137.09°, 
and Md = 289.23°. For each synthetic measurement (vobs), 
we compute the RV value at that time (vmod), then add two 
noise terms: one drawn from a standard normal scaled by 
the measurement uncertainty aobs and a similar term but 
scaled by the jitter value corresponding of that observation, 
i.e. Vobs — n^nod “b A7(0, 1) X CJohs “b A7(0, 1) X CjH- 

These datasets were then used as simulated input obser¬ 
vations for RUN DMC, and a full differential evolution MCMG 
was performed on each synthetic data set (in the same man¬ 
ner as was performed on the real data, and described in 
m\ and [4]). We performed three realizations of the RV time 
series for both Schaos-.so and Siieg-.20- 

A2 Results for Synthetic Planetary Systems 

In right panel of Figure [M\ we show the (pLapiace libra¬ 
tion amplitude distribution for the three realizations of the 
RV time series based on Schaoa-.so- The chaotic nature of 
Schaos ;8o is shown in Figure HD left where the periodicity 
and peak-to-peak variation are not constant in time. There 
is a rough variation of ~160° corresponding to an amplitude 
of ~80°. For each realization, we find libration amplitudes 
of 97±ig (blue dash-dotted), 96±2i (red dashed), and 91±i7 
(green solid) degrees based on 1,000 solutions each. Only a 
few solutions from each of these realizations were unstable 
over the short integration baseline (1 kyr). These estimates 
are qualitatively consistent with the input amplitude. 

In right panel of Figure IA2I we show the 4>Lapiace li¬ 
bration amplitude distribution for the three realizations of 
the RV time series based on Schaos-.so- A significant fraction 
of our posterior sample for each realization was not dynam¬ 
ically stable. Specifically, only 548, 496, and 673 systems 
remained after the 1 kyr integration. We find that the RV 
analysis of the synthetic data overestimated the outer-most 
planet’s mass by ~20%, which could be reason for such a 
large fraction of unstable systems. We find libration ampli¬ 
tudes of 25.6±5;3 (blue dash-dotted), 22.9±5;5, (red dashed) 
and 19.5±3;7 (green solid) degrees based on the remaining 
short-term stable systems. Again, these values are qualita¬ 
tively consistent with the input amplitude of ~20°. 

It Is clear from these results that our RUN DMC algorithm 
does not significantly bias the recovered libration amplitudes 
for the 4:2:1 resonance of the outer three planets in the GJ 
876 system. Thus, we believe that our best-fit libration am¬ 
plitude for the real data (® is accurate and the GJ 876 sys¬ 
tem is in a chaotically librating state, with a high-amplitude 
libration. 

A3 Estimation of Lyapunov times for GJ 876 
Orbital Solutions 

iBatygin et al.l (l2015l l suggests that the GJ 876 system is 
chaotic, with a characteristic timescale (or Lyapunov time, 
defined below) for the chaos of roughly 14 years. This the¬ 
oretical study ignored the innermost planet and treated 
the system as coplanar. An estimate of the Lyapunov time 
for the 4-planet, coplanar radial velocity orbital solution of 
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Figur e Al. Radial velocity analysis of a synthetic system generated from a turbulent migration model, 3(^/1,0,05:80 5 from iBatvgin et alJ 
1120151 ). Left: The chaotic evolution of (pLaplace over the course of 100 years for model Schaos-.SO- R-ight: The distribution of the (pLaplace 
libration amplitude based on three realizations (green solid, red dashed, blue dash-dotted) of the RV data generated from the synthetic 
system, Schaos-.SO- 



Time [kyr] 



Figure A2. Radial velocity analysis of a synthetic system generated from a smooth migration model, S/{eg. 2 o. Left: The evolution of 
the absolute value of the cj>Laplace over the course of 4 X 10® years for model Si{e 9 : 20 - R-ight: The distribution of the ipLaplace libration 
amplitude based on three realizations (green solid, red dashed, blue dash-dotted) of the RV data generated from the synthetic system, 
SReg:20- 


iR.ivera et al.l ll20inl) was c onsistent with the analytic esti¬ 
mate. Barnes et al.l ll2015l) suggests that mutually inclined 
systems in or near a MMR can also exhibit chaotic evolution. 

Here we study the nature of three sets of solutions: 
first, the 1,000 long-term stable, 3-d orbital models de¬ 
scribed in SjS] 1,000 coplanar models described in ||4l and 
1,000 coplanar models based on the SReg-. 20 - We deter¬ 
mine whether or not the orbits are chaotic by evolving si¬ 
multaneously the standard gravitational equations of mo¬ 
tion and the variational equations of motion, which yields 
an estimate of the Lyap u nov t ime of a trajectory (e.g. 
iLichtenberg fc LiebermanI (Il992l )). The variational equa¬ 
tions govern the behavior of small perturbations to an orbit 
and therefore can be used to study how perturbations evolve 
in time. For chaotic orbits, small perturbations of length D 


grow exponentially as D ~ with a characteristic time 
r, which in the limit as f —>■ oo is defined as the (mini¬ 
mum) Lyapunov time. Our finite time integrations are used 
to estimate this Lyapunov time as f/inai/log [(D(t/inai)], 
where tfinal is the total integration time, D{t = 0) = 1, and 
D{tfinal) is the total length of the “perturbation” at the 
end of the integration. Note that D in principle can become 
very large, but since the variational equations are linear in 
the components of D, the absolute length of D need not be 
small for the variational equations to apply. If an orbit is 
regular, the reported Lyapunov time will be comparable to 
the integration time, though integrations cannot prove an 
orbit is regular. These integrations must therefore be car¬ 
ried out for long enough such that the chaotic and regular 
orbits have markedly different reported Lyapunov times. 
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Figure A3. The distribution of estimated Lyapunov times for 
the GJ-876 system. These estimates are the result of integrations 
of the set of 10^ 3-d orbital solutions fit to the radial velocity 
data, employing a time step of 0.14 days (red) and with a time 
step of 0.014 days (black), and for the 1,000 solutions fit to the 
data assuming coplanar orbits, which employed a time step of 
0.14 days (blue). 


We employed a Wisdom-Holman mapping in canonical 
astrocentric coordinates to integrate b oth the equations of 
motio n and the variational equations l|Wisdom fc HolmarJ 
Il99lh . A third-order symplectic corrector wa s implemented 
to im p rove the accuracy of these integrat ions (I Wisdom et al.l 
1 19961 : IChambei^ 1 19991 : IWisdomI 120061 ). We used a sim¬ 
ple prescription for the effects of general relativity which 
leads to precession of t he orbits with the correct timescale 
llMilani fc Nobilil 1 19831 ). As we will show, the Lyapunov 
times of the orbits were significantly shorter than the 
timescale of precession due to general relativity, and so this 
approximate version of general relativity is sufficient for a 
good estimate of the Lyapunov times for the orbits studied. 

The Wisdom-Holman integrator can be used with a 
time step as large as a tenth or twentieth of the shortest 
orbital timescale in th e system llWisdom fc HolmarJ Il99^ : 
iRauch fc Holmanlll999l ). For GJ 876, this corresponds to the 
time needed to resolve the pericenter passage of the inner¬ 
most planet. We estimate this using the orbital period of the 
innermost planet as if its semimajor axis was equal to the 
pericenter distance, Puff ^ Porb{l — ~ 1.6 days, where 

e « 0.1 and Porb ~ 2 days. We used a time step of either 
0.14 or 0.014 days when integrating the set of 3-d solutions, 
and, as discussed below, these yield similar estimates of the 
Lyapunov times. For both coplanar set of solutions, we used 
a time step of 0.14 days. The maximum fractional energy 
error /S.E/E was ~ 10~® for the 3-d set of orbits with a 
time step of 0.14 days, ~ 10“^° for the same set with a time 
step of 0.014 days, and ~ 10“^ for the coplanar set of orbits. 
These integrations lasted 10^ years. 

In Figure IA3I we show the distribution of Lyapunov 
times for the 3-d set of solutions resulting from the inte¬ 
grations employing both timesteps, and the distribution of 
Lyapunov times for the coplanar set of solutions. These in¬ 


tegrations all agree that the motion of these orbits is chaotic 
with a Lyapunov time of roughly 10 ye a rs, con sistent with 
the analytic estimate of iBatvein et al] (l2015h . The small 
peak at 10® years in the distribution of Lyapunov times for 
the coplanar orbital solutions to the data corresponds to or¬ 
bits which might be regular. However, the vast majority of 
coplanar orbits studied were chaotic. The close agreement 
between Lyapunov times estimated for the 3-d and copla¬ 
nar se ts is also consistent with the analysis of lBatvein et al.l 
l|2015l) and iMartf et al.l (l2013l ) in that the chaotic motion is 
captured by a coplanar approximation of the true system. 

During the integrations with a time step of 0.14 days, 
~10 orbits in each set (3-d vs. coplanar) were flagged as un¬ 
stable. With the smaller time step of 0.014 days, no orbit in 
the set of 3-d orbits was flagged as unstable. This suggests 
that the smaller time step might be necessary for study¬ 
ing the long-term stability of these orbits using a Wisdom- 
Holman integrator. However, since the two time steps agree 
on the distribution of Lyapunov times, we believe these re¬ 
sults are robust. Lastly, we do not show the results of the 
coplanar set of solutions fit to synthetic data, since our inte¬ 
grations suggested that nearly all of these orbits where not 
only chaotic but also showed instability on short (10^ year) 
timescales. We defer any further analysis of the stability of 
these synthetic solutions to future work. 
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